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Abstract 

We aimed to use finite element method to simulate the unique behaviors of liquid crystal 
elastomer, such as semi-soft elasticity, stripe domain instabilities etc. We started from an energy 
functional with the 2D Bladon- Warner- Terentjev stored energy of elastomer, the Oseen-Frank 
energy of liquid crystals, plus the penalty terms for the incompressibility constraint on the dis- 
placement, and the unity constraint on the director. Then we applied variational principles to 
get the differential equations. Next we used mixed finite element method to do the numerical 
simulation. The existence, uniqueness, well-posedness and convergence of the numerical methods 
were investigated. The semi-soft elasticity was observed, and can be related to the rotation of 
the directors. The stripe domain phenomenon, however, wasn't observed. This might due to the 
relative coarse mesh we have used. 

1 Introduction 

Nematic liquid crystal elastomer is a relatively new kind of elastic material which combines proper- 
ties from both incompressible elasticity and rod-like liquid crystals. It has some unique behaviors 
such as the stripe-domain phenomenon and the semi-soft elasticity. These phenomena are usually 
observed in the "clamped-pulling" experiement, in which a piece of rectangular shape liquid crys- 
tal elastomer is clamped at two ends and pulled along the direction perpendicular to the initial 
alignment of the directors. 

Mitchell et al. |17j did the "clamped-pulling" experiement for acrylate-based monodomain 
networks. They found that the directors rotate in a discontinuous way: after a critical strain, 
the directors switch to the direction perpendicular to the original alignment. On the other hand, 
Kundler and Finkelmann [13] did the same experiment for polysiloxane liquid crystal elastomer. 
They found that the directors rotate in a continuous manner. They also observed the interesting 
stripe domain phenomenon. Later, Zubarev and Finkelmann |20j re-investigated the "clamped- 
pulling" experiment for polysiloxane liquid crystal elastomers, to investigate the role of aspect 
ratio of the rectangular shape monodomain in the formation of stripe domains. Figure [T] shows the 
process of the "clamped-pulling" experiement in the case that the aspect ratio is AR = 1. In this 
graph, the pulling direction is vertical, while the initial alignment of the directors is horizontal. 
The size of the elastomer was 5 x 5mm 2 , with thickness 50 — 300 ± 5 /im. They found that at 
some critical elongation factor A = 1.1, an opaque region started to emerge in the center of the 
elastomer. Using X-rays, they found that this region was made of stripes of about 15/im in width. 
In each stripe, the directors were aligned uniformly, while across the stripes, the directors were 
aligned in a zig-zag way, as shown in the second picture of Figure [l] When the elastomer was 
pulled further, the opaque region was broken into two regions symmetric about the center (the 



Figure 1: The stripe domain phenomenon. Reproduced from |20j 



third picture of Figure [T]) . And these two regions moved closer to the two ends as the elastomer 
was pulled further. Eventually at A = 1.3 the two stripe domain regions reached the two ends 
(the last picture of Figure [I]) and in this last stage, most of the directors at the center had rotated 
90 degrees, and were now aligned in the pulling direction (vertical direction). 

Another interesting phenomenon of liquid crystal elastomers is the so-called "semi-soft elas- 
ticity". Figure [2] shows the stress-strain graph of three different kinds of liquid crystal elastomers 
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Figure 2: The semi-softness of liquid crystal elastomer. Reproduced from [19J. 

in the "clamped-pulling" experiment (P3J [7]). We can see that at the beginning, the nominal 
stress increases linearly with the strain, just like ordinary elastic materials. However, after the 
strain reaches a critical point, the stress-strain curve becomes relatively "flat", which we call the 
"plateau region" . When the strain falls in this region, relatively small stress can induce relatively 
large deformation, which means the elastomer is "soft". We can also see from Figure [2] that when 
the strain continues to increase and reaches another critical value, the stress-strain curve becomes 
linear again. 

The stored energy of liquid crystal elastomer proposed by Bladon, Warner and Terentjev 
(BTW) is [2] 

W(F, n)=n (|F| 2 - (1 - a)\F T n\ 2 - 3a 1/3 ) (1) 

where F = I + Vtt(X) is the deformation gradient of the elastomer network, while n = n(X) is 
a unit vector denoting the orientation of the rod-like liquid crystal polymers, and X denotes an 
arbitrary point in the reference configuration. We assume that the elastomer is incompressible, 
which is equivalent to det(F) = 1. Also fx is the elasticity constant, and a € (0, 1) is a constant 
measuring the significance of the interaction between the displacement u and the orientation n. 
Notice that in the limit a = 1, the BTW energy degenerates to the stored energy of incompressible 
nco-Hookean materials. The closer the value a is to 1, the smaller the interaction between the 
displacement u and the orientation n. Similarly in 2D, the BTW energy can be defined as 

W(F, n)=n (\F\ 2 - (1 - a)\F T n\ 2 - 2a 1/2 ) (2) 
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The only difference is that the constant term is changed from 3a 1 / 3 to 2a 1 / 2 . 

Proposition 1. The BTW energy W(F,n) in equation |7p and are always non-negative. 

• In 3D, W(F,n) = if and only if eig(FF T ) = {a 1 / 3 , a 1 / 3 , a~ 2 / 3 } and n is an eigenvector 
of the eigenvalue a~ 2 / 3 . 

• In 2D, W(F,n) = if and only if eig(FF T ) = {a 1 / 2 , a -1 / 2 } and n is an eigenvector of the 
eigenvalue a" 1 / 2 . 

Proposition 2. • In 3D, assume \{t) = \{t) = a 1 / 6 (l - t) + a' 1/3 t with < t < 1. When 
there is no shearing and 

/a 1 / 6 \ 

F = A(t) (3) 

\ a- l ^/\{t)J 

the BTW energy |7J) is zero /or t — 0,1 and strictly positive for < t < 1. However, when 
there is shearing and 

fa 1 / 6 \ 

F = I A(t) ±5 (4) 
V a-W/X®) 

with 

5 = ^/aV 3 + a- 2 / 3 - (A 2 + a- 1 / 3 /* 2 ), (5) 

ine BTW energy |7p can 6e zero /or < t < 1. 

• In assume X(t) = a 1 / 4 (l — i) + a _1 / 4 i iwt/j < t < 1. When there is no shearing and 

i/ie BTW energy is zero /or £ = 0, 1 and is strictly positive for < i < 1. However, 
when there is shearing and 



V V-M*) 

with 

S = ^/a 1 /2 + a -i/2_( A 2 + A-2), (8) 

ine BTW energy |IJ) can 6e zero /or < t < 1. 
Corollary 3. //we iafce 

= min W B ™(n,F), (9) 
|n|=l 



V (A? + A! + aA 2 - 3a 1 / 3 ) in 3D 
fi (A 2 + aXj - 2a 1 / 2 ) in 2D 



(10) 



inen W(.F) is a non-convex function of F. 

Proof. Take < i < 1, and X(t) and F± as defined Proposition^ Then we have 

W(0.5F + + 0.5F_) > 

while 

0.5W(F+) + 0.5VF(F_) = 
Therefore W(F) is a non-convex function of F. □ 
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Proposition [2] tells us that introducing shearing might lower the BTW energy or @. On 
the other hand, because of the constraint of the clamps, global shearing is not possible, therefore, 
local shearing is developed in a zig-zag way, and that's why stripe domains occur. 

In [5], DeSimone et al. did the finite element study of the clamped pulling of liquid crystal 
elastomer. They started from the 3D BTW energy , and got rid of the variable n by taking 
it to be always n = n(F), the eigenvector of the largest eigenvalue a -2 / 3 of the matrix FF T . So 
they got 

W(F) = n (\\{F) + Xl(F) + a\l(F) - 3a 1 / 3 ) (11) 



It turned out that the energy function (11) is no longer a convex function of F (Corollary pj). 



Therefore, they replaced it by its quasi-convex envelope 

W qc (F) := inf i i^T / W(Vy(x))dx : = Fr on Oil } (1.2) 



dct(Vj/(x)) = l 
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They pointed out, that the use of W qc (F) in the numerical computations allows one to resolve only 
the macroscopic length scale, with the (possibly infinitesimal) microscopic scale already accounted 
for in W qc (F). 

It turned out that their approach was very successful, in that they successfully captured the 
emerging and migrating of the stripe domain regions, in exactly the same way as those observed 
experimentally in [20) . However, their approach is limited for several reasons: 

• Although they were able to tell which part of the region is a "stripe-domain region" , they 
couldn't tell how many stripes lie in that region. Actually, their model permits infinitely 
many stripes in that region, which is unphysical. 

• They were only able to observe the "soft elasticity" , not the "semi-soft elasticity" . That is, 
the plateau region occurs immediately upon the pulling (Figure [3]), instead of after some 
critical strain as observed experimentally. 



1 1.2 1.4 1.6 




1 1.2 1.4 1.6 l.B 2 



strain 

Figure 3: The stress-strain curve from Desimone's numerical simulation. Soft elasticity instead of 
semi-soft elasticity was observed. Reproduced from jS]. 

In this work, we aim to resolve the above issues. We combine the 2D BTW model ^ and the 
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Oseen-Frank model ([H]), and get the following energy functional for liquid crystal elastomers: 
H(u,n)= / c(| J F| 2 -(l-a)|F T n| 2 )+6|Vn| 2 



f ■ u — g ■ uda 

n Jr 2 

where / is the body force, while g is the traction force on the T 2 part of the boundary dft. We 
assume u £ and satisfies det(7 + Vit) = 1 a.e. in fi; while n £ H 1 ^), and satisfies 

|n| = 1 a.e. in fl. Let L be the characteristic length scale of the domain ft, we get the following 
non-dimensionalized energy 



n(u,n) = f (\F\ 2 - (1 - a)\F T n\ 2 ) +b\\7n\ 
Jn 



f ■ u — / g ■ uda 

where ft is congruent with ft with characteristic length scale 1, and T2 C <9ft is the image of T2. 
Also6= o | y ,/ = i,a n dg=a. 

The rest of the paper is organized as the following. In section [2j we investigated the properties 
of the continuous problem, such as the existence of minimizer, the equilibrium equation and its 
linearization, the stress-free state, the well-posedness of the linearized equation etc. In section 
[3j we investigated problems related to the discretization, such as the existence of minimizer, 
the equilibrium equation and its linearization, the well-posedness of the linearized equation, and 
the existence and uniqueness of the Lagrange multipliers etc. In section |4j we presented the 
simulation results for the "clamped-pulling" experiment, including the results of the inf-sup tests, 
the convergence rates, and the stress-strain curve etc. Finally, in section [5] we discussed future 
directions. 

2 The Continuous Problem 

2.1 Existence of minimizer 

Let 

£T 1 (f7,S' 1 ) = {n £ H 1 ^,^ 2 ) : \n\ = 1 almost everywhere in ft} (13) 

and 

K = {u£ H l (Q,, K 2 ) : det(J + Vit) = 1 almost everywhere in ft} (14) 
Define the admissible set 

A(u , n ) = {u £ K,n £ -ff 1 (ft, S 1 ) : u = u on r , n = n on Ti} (15) 

Define the energy functional 

n(u,n)= f (|F| 2 -(l-a)|F T n| 2 )+6|Vn| 2 
Jn 



f u 9 uda (16) 

Jn Jr 2 

where F — I + Vii. Then our problem is 

Find (u,n) £ A(uo,no) minimizing n(w,n). (17) 
Lemma 4. Assume \n\ — 1, we have 

(|F| 2 -(l-a)|F T n| 2 )>a|F| 2 (18) 
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Proof. Let X\ < \\ be the two eigenvalues of the real symmetric matrix FF T . Since n is in 
H 1 ^^ 1 ), by Proposition [lj wc have 



{\F\ 2 -(l-a)\Fn\ 2 )>\\ + a\l 

2 i \2\ 



>«(Af + A^) 
= atr(FF T ) 
= a|F| 2 

□ 

Lemma 5. Assume \n\ = 1, and < a < 1. T/ie function 

L(F,n) = \F\ 2 -(l-a)\F T n\ 2 (19) 

is convex in F. 
Proof. Let 



Then we have 



A(n) = I - {1 - a)nn T (20) 
L = tr(FAF T ) (21) 



Then for any Fi,F 2 £ Mr, and < a < 1, we have 

[aL(Fi) + (1 - a)L(F 2 )] - L(aFi + (1 - a)-F a ) 
= crtr(FiAF\ T ) + (1 - a)tr(F 2 AF 2 T ) 

- tr[(aF a + (1 - a)F 2 )A{aF 1 + (1 - a)F 2 ) T ] 
= a(l - o)tr[(fi - F 2 )A(F X - F 2 )] 
= a(l - a) (\F t - F 2 | 2 - (1 - o)|(fi - F 2 ) T n| 2 ) 

> 2a(l - a)V^ 

> 

where we have used Proposition [T] Therefore L is a convex function of F. □ 
Remark It's easy to see that Lemma [4] and Lemma [5] holds for 3D, too. 

Theorem 6. Let 51 be a nonempty, bounded, open subset ofR d . 

• If d = 2, suppose we have Uk u in W 1 ^ with s > | ; then we have det(J + Vitfe) — > 
det(/ + Vu) in V'(Vt). 

• J/d = 3, 

— suppose we have Uk — 1 u in W 1,s with s > |, £/ien we /iave adj(I + Vttfe)^ — > adj(J + 

— suppose we have Uk — 1 m in W^ 1 ^, and adj(/ + Vuj) — 1 adj(7 + Vii) in L 9 (57;Mt 3 ) wii/i 
s > 1, q > 1 and ~ + A < | } t/ien we have det(7 + Vw fe ) -> det(J + Vu) in X>'(fi). 

For the proof, please see John Ball pQ. 
Theorem 7. There exists solution to the problem Xl 
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> 



Proof. Let m be the infimum of II on A(uo,no), and let (life, rife) £ A(uq, no) be a minimizing 
sequence of II. Obviously m < +00. Thus H(u k , Tik) is bounded above. By Lemma|4j we have 

Tl(u k ,n k )> / a\(I + Vu k )\ 2 + b\Vn k \ 2 dx 

~ ll/llL 2 (0)ll M fc||L 2 (J2) ~ llfl'llL 2 (r 2 )ll'"A ; !|L2( r2 ) 

/ a\(I + Wu k )\ 2 + b\Wn k \ 2 dx 
Jn 

\\f\\h(n) +4uk\\h ( n^j - (^\\9\\h(r 2 )+4uk\\ 2 L 2 (r2 
> [ d\(I + Vu k )\ 2 + C 2 \Vn k \ 2 dx - C 3 

where e > is small and Cj > 0, i = 1,2, 3 are constants and we have applied the generalized 
Poincare inequality ({4], p281) and the Trace Theorem ([9], p258) in the last step. Therefore we 
have F k = I + Vu k and Vn k are bounded in L 2 . Now we have that V(u k — u ) is bounded in 
L 2 , by the Poincare inequality, we have u k is bounded in H 1 . On the other hand, since Vrij. is 
bounded in L? and n k is in H 1 ^!, S 1 ), we can get that n k is bounded in H 1 . Now since H 1 is a 
reflexive Banach space and u k and n k are bounded in H 1 , we can find a subsequence of u k and a 
subsequence of n k such that they are weakly convergent in H . We still denote them as (u k , n k ), 
and assume u k — v it, ra^ — 1 n. 

Now that itfc — ^ u in i? 1 , by Theorem [6j we have det(I + Vu k ) — » det(I + Vm) in ©'(fi). 
Since det(7 + Vufc) = 1 a.e., we have det(I + Vu) = 1 a.e. in £1 Q On the other hand, weak 
convergence in H 1 implies strong convergence in L 2 [^] thus we can find a subsequence of n k that 
converges point-wise almost everywhere. Therefore we have \n\ = 1 almost everywhere. Finally 
since u k — Uq e -^oirv an< ^ ^o|r ^ s a c l° se d linear subspace of H , by the Mazur's Theorem, it's 
weakly closed. Therefore it — «o is also in Hh^ , which means u — Uq on Tq. Similarly we have 
n = no on Y\. Therefore we have (u,n) e A(uo,n,Q). 

By Lemma [5j the following function 

L(F, n, P) = {\F\ 2 - (1 - a)\F T n\ 2 ) + b\P\ 2 

is a convex function of F and P. Therefore by Theorem 1 of section 8.2 of jS], II is weakly lower 
semi-continuous. Thus we have 

n(ii, n) < liminf n(ttfe, rife) 
= m 

Since m is the infimum of II on A, we conclude that 

n(u,n) = m (22) 
That is, (it, n) is the minimizer of LI on A. □ 

Remark Theorem[6]is crucial in this existence proof. Since in the 3D case, we need the additional 
condition that a.dj(F k ) is weakly convergent to get the convergence of det(F k ), the above proof 
cannot be directly extended to 3D. 



1 This is because, by definition, we have 

(det(J + Viife), 4>) -» (det(J + Vu), <j>) 
in M for any <f> £ T>(Q). Since det(J + Vu k ) — 1 a.e. in f2 for any k, we have 

(det(J + Vm) - 1, </>} = 0, V(fe»(Sl) 

Thus det(7 + Vit) = 1 a.e. in Q. 

2 This is because the embedding / : W 1 ^ — > L p is compact for 1 < p < 00 ( 9 , p274), while for any compact operator 
A : V —tW with V and W Banach spaces, u k — =■ u in V implies Au k — > Au in W ([4], Theorem 7.1-5 on p348). 
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2.2 Equilibrium equation and its linearized system 

Constrained minimization problems can usually be reduced to unconstrained minimization prob- 
lems by introducing Lagrange multipliers. We consider the following non-dimensionalized energy 
functional 

£(u,p, n, A) = / (\F\ 2 - (1 - a)|F T n| 2 ) + 6|Vn| 2 

-p{dctF-l) + \{\n\ 2 -1)- [ f -u- [ gu (23) 

Jn Jr 2 

where p is the Lagrange multiplier for the incompressibility constraint det(J+Vu) = 1, which can 
be interpreted as "pressure", while A is the Lagrange multiplier for the unity constraint \n\ = 1. 

Assume (it, n,p, A) minimizes the non-dimensionalized energy (23). Then for any test function 
v, the function £ (e) = £{u + ev,p, n, A) has a minimum at e = 0. Thus we have 

0-|(0) (24) 
from which we can get the following equation 

0= / 2 (F : Vv - (1 - a){F T n, Vv T n)) - p^-^ ■ Vv 
Jn ob 

f v- g ■ vda 



By similarly taking the variations n — »• n + sm, p — > p + sq, or A — > X + sp,, we get the following 
set of Euler-Lagrange equations 

0= / 2 (F : Vv - (1 - a)(F T n, Vv T n)) - p^^ : Vv 
Jn 

- fv- g ■ vda 

Jn Jr 2 

0= / -2(1 - a)(F T n,F T m) + 26(Vm, Vn) + 2A(n,m) 
Jn 

= / -g(detF-l) 
0= / /*((», n)-l) 

in which we seek solutions (u,n,p, A) € Hgi ro x H^,i ri x L 2 (£l) x L 2 (f2). Correspondingly, the 
test functions (v,m,q, p) are in x Hg, ri x L 2 (Vl) x L 2 (il). Here -ffg|p is defined as 

ff fl Vo = { ,;eifl ( Q ; R )^lro=5} 

and H^|p o is the corresponding vector version. 

We linearize the Euler-Lagrange equation around a solution (u, n,p, A), and get the following 
linearized equations: 

a 1 (w,v)+a 2 (l,v)+bi(o,v)=Li(v) Vt> € V (25) 

a 2 {m,w) +a 3 (l,m) + b 2 (7, m) = L 2 (m) VmeM (26) 

6i(g,«) = L 3 (g) V<? e P (27) 

b 2 (fx,l) = L A (fi) VpeA (28) 
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where V = Hj, ro , M = Hj, ri , P = L 2 (fl), A = H^, and (to, 1, 0,7) eVxMxPxAis the 
change in the solution, (v, m, q,/i) eVxMxPxAisa test function. Here H^^ is the dual space 
of Hh Ti . We'll abuse the notation and use || • ||_i to also denote the norm in the Hf* space. The 
bilinear forms in the linearized system are defined in the following equations 

ai(w,v)= / 2(Vw, Vv) - 2(1 - a)(Vw T n, Vv T n) 
Jn 

,d 2 det . „ ,„„, 
-K-^-Vt£>):Vo (29) 

a 2 (m, v)= I -2(1 - a){F T m, V« T n) - 2(1 - a){F T n, V« T m) (30) 
Jn 

a 3 (l,m)= I -2(1 -a)(F T l,F T m) + 26(Vrn,VZ) + 2\(l,m) (31) 
Jn 

bx(q,w) = / ~q—:Vw (32) 



b 2 ( P J) = [ 2/i(i,n) (33) 
Jn 



2.3 The stress-free state 



From the Euler-Lagrange equations, we can get the following strong form partial differential 
equations: 

-diva = / in fl (34) 

fediv(Vn) + (1 - a)n T FF T - Xn T = in fl (35) 

det F = 1 in (36) 

(n, n) = l in ft (37) 

and the following natural boundary conditions: 

<jv = g on <9£!\ro 
Vnv = on <9Q\:Ti 

where the Piola-Kirchhoff stress tensor is 

a = 2(l-(l-a)nn T )F-p^ (38) 

Notice that if we choose the displacement to be ti e and the director to be n = (0, 1) T , then 
we have 

f 2-p 
2a - p. 



If there is zero body force, equation (34) implies that p must be a constant. Then there is no way 
to make this a zero if a 7^ 1. However, it can be verified that the stress-free state can be achieved 
with 

'a 1 '* 



and 

n=(0,l) T 

p = 2\fa 

A = (l-o)/Vo 
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Remark The above observation implies that the reference configuration (strain-free state) for 
the BTW model is different from the stress-free configuration. It can be verified that if we take 
the reference state to be the stress-free state instead, then the BTW model becomes the so-called 
"general BTW" model 

W B tw = ~ f J,tv(L F T L- 1 F) (40) 

where 

L(n) = a-^ 2 nn T + a^ 2 (I - nn T ) 
is the so-called step tensor, and Lq = L(tiq) is the step tensor at the stress-free state. 



2.4 Well-posedness of the linearized system 



By adding ([25]) and ( |26| together, and adding ( 27 )-( 28 ) together, we can reduce the linearized 
system into a standard saddle point system 



a(w,v) + b(o,v) = Lx(v) VieY 
b(q,w)=L 2 (q) VqeP 

where V = ¥ x M, P = P x A, and 

a(w, v) = a±(w, v) + (12(1, v) + 02(171, v) + 0,3(1, m) 
b(q,v) = bi(q,v) + b 2 ((J,,m) 



(41) 
(42) 



(43) 
(44) 



The well-posedness of saddle point systems are well established. Here we quote the so-called 
Ladyzenskaya-Babuska-Brezzi theorem from [IS] , 

Theorem 8 (Ladyzenskaya-Babuska-Brezzi). Consider the following saddle point problem 

a(u,v) + b(p,v) = L v (v) WeV,u€V (45) 
b(q, u) = L P (q) V 9 € ¥,p € P (46) 

with V and P given Hilbert spaces, Ly and Lp belonging to V and P' respectively. Moreover, a 
and b are continuous bilinear forms defined on V x V and P x V, respectively. Define the operators 



v i-> 38v such that (SBv, q) = b(q, v) V?eP 

w <— > srfw such that (srfw, v) = a(w, v) Vu G Ker^ 

Then the operator S3 is onto if and only if the spaces V and P satisfy the following inf-sup 
condition: 

inf sup b(q,v)>/3>0 (47) 

geP,||g||=i „<=v,|M|=i 

Moreover, the mixed problem is well-posed if and only if 3§ is onto and £/ is invertible. 
Remark The operator srf is invertible if and only if the following inf-sup condition is satisfied 

a(w, v) 



inf sup 



> a > 



w v 



For some mixed system, we can prove the so called "ellipticity" condition 

a(v, v) 



inf 

uEKcrJ \\V\\ 



> a > 0, 



(48) 



(49) 



which is a stronger condition than the inf-sup condition ( 48 ) . Thus the ellipticity condition ( 49 ) 



system (45)- (46 1 



together with (47 1, will give us a sufficient condition for the well-posedness of the saddle point 
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Theorem 9. The inf-sup condition for b{q, v) = b\{q, v) + b 2 ([j, 1 tn) is satisfied if and only if the 
corresponding inf-sup conditions for bi(q,v) and b 2 {/i,m) are satisfied. 

Proof. First if the inf-sup condition for b{q, v) is satisfied, then by Theorem [8j we know that the 
operator 

I:Vxi4P'xA' 

(v, m) 1 — ^ £${v, m) such that {38 {v, m), (q, /i)) = b 1 {q, v) + b 2 {n, m) V(g, fi) e P x A 
is onto. Therefore, the operator 
38 x :V P' 

v (— > such that {38iV,q) = h(q,v) Vg e P 

and the operator 

^2 :M -> A' 

m 1 ^ 38 2 m such that {382m, A*) — ^(M; wi) V/ie A 

are both onto. Thus by Theorem [8j the inf-sup conditions for 61 (g, u) and b 2 {/j,,m) are satisfied. 

Conversely, if the inf-sup conditions bi{q,v) and b 2 (/j,,m) are satisfied, then the operators ^1 
and ^2 are both onto. Thus the operator 38 is also onto. Therefore we have the inf-sup condition 
for b(q, v) — h{q,v) + b 2 {^ 1 m). □ 

Therefore, to verify the inf-sup condition for b{q, v) = b\ (q, v) + b 2 (fi, m) , it's enough to verify 
the inf-sup conditions for bi{q, v) and b 2 (/i, m) individually. Notice that the bilinear form b\(q, v) 
in (32) is exactly the one in the incompressible elasticity [TB], while 62 (a*, in (33) is exactly 
the one in the harmonic map problem \V2\ . 

For the inf-sup condition for b±{q, v), it's well-known that it's satisfied at the strain-free state, 
where F = I, and the inf-sup condition is reduced to the one in the Stokes problem: 

inf sup #S^>/3r>0 (50) 
96L 2 (n)„ eH i m ||g|| |M|i 

Since the stress-free state has constant F matrix, by change of variables, it's easy to verify that 
the inf-sup condition for b\(q, v) is satisfied at the stress-free state, as well. In the general case 
that u ^ and F is not a constant, analytical verification of the inf-sup condition for b\ {q, v) can 
be very difficult. 

As for the inf-sup condition for b 2 (/i, m), a slight modification of the proof in [TJ] gives us the 
following theorem: 

Theorem 10. Assume n e Hj, ri (f2) |~| W 1,oa {fl), then the second inf-sup condition for b 2 {fi, m) 
is satisfied, that is 

inf sup ,{ 2 "' ™!. m) > [3 2 > 0. (51) 

peH- i 1 (o) TneH i |ri (o) l|n»l|i||Mll-i 

Finally, for the ellipticity condition for the bilinear form a{w, v), it's generally very complicated 
to verify due to the complexity of the expressions of a\{-, •), a 2 {-, ■) and a%{-, •). However, it can 
be verified that at the stress-free state, the ellipticity condition for a(w,v) is not satisfied. This 
doesn't necessarily mean that the linearized system ( 25 )-( 28 ) is not well-posed, though. After all, 
as remarked before, the ellipticity condition is a sufficient condition, not a necessary condition. 
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3 Discretization 



In the elastomer problem, we have the following variables to solve: 

• The displacement vector field u and the pressure p, which is also the Lagrange multiplier 
for the incompressibility constraint det(F) — 1 = 0, 

• The director vector field n and the Lagrange multiplier A for the unity constraint |n| 2 — 1 = 0. 

The first pair (u,p) is similar to those in the incompressible elasticity, and we'll use the Taylor- 
Hood element P2 — Pi for (u,p). That is continuous piecewise quadratic finite element for it, and 
continuous piecewise linear finite element for p. The Taylor- Hoold element has been proved (|18j) 
to be stable at least at the strain-free state u = 0. The second pair (n, A) is similar to those in 
the harmonic map problem, and we will use piecewise linear finite element for both n and A, as 
Winther et al. did in [T5]. Notice that it's crucial to impose Dirichlet boundary conditions for n 
and A at the same boundary to make sure that |n| = 1 is satisfied at all the mesh nodes. 

Let Vh denote the space of continuous piecewise linear functions and Vh, g \r = { v Vh H H 1 : 
v = g on r }. Let V/j and ~Vh,g\r be the corresponding vector version. Let be the nodal 
interpolation operators onto the spaces Vh and Vfe. 

Let Wh denote the space of continuous piecewise quadratic functions and W^gWo = i w 
Wh H H 1 : w = g on Lq}. Let and W; l g |r be the corresponding vector version. 



3.1 Existence of the discrete minimization problem 

Let 

K h = {u h g W h , olra + u oh , / q h (det(I + Vu h ) - l)dx = 0,Vq h G V h } (52) 



N h = {n h G V M | ri +n oh , / fj, h ^h{\n h \ - l)dx = 0, V/ift G Vh,o|rJ (53) 

Jn 

Define the admissible set 

A(u oh , n oh ) = K h x N h (54) 

Notice that rih € Nh if and only if the function 7T/i(|to/j| 2 — 1) G Vi l: o\r 1 is identically 0, which 
means \rih\ = 1 at all the mesh nodes. 

Let <pj,j = !,••■ , N be a basis of Vh, and = !,••• , M be a basis of 14, |ri- And define 



g 3 [Uh , n h ) - i ^_^ ( | n/i |2 _ 1)da . AT + 1 < j < iV + M 



Then gj is a continuous function on (W^Qiro +U0/1) x (V^ t o|ri + n oh)- Therefore ^(M ^,n,o^) 
can be written as the intersection of reciprocal images of of the continuous functions gj . Thus 
it's closed in (W M |r + u oh) x (V fe|0 | ri +n oh )- 
Define the energy functional 

n(u,n)= / (|F| 2 -(l-a)|F T n| 2 )+6|Vn| 2 
Jn 

f ■ u g ■ uda (56) 



where F = I + Vm. Then the discrete formulation of the minimization problem is 

Find (v,h,Tih) G A(uoh,Tioh), minimizing II on A(u oh , n oh ). (57) 
Lemma 11. Assume n £ Nh and < a < 1, then for any matrix F G M 2x2 ; we have 

(\F\ 2 ~(l-a)\F T n\ 2 )>a\F\ 2 (58) 
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Proof. Take any point x £ fl, suppose it's in the triangle APiPapj. Since n £ Nh, we have 

n(x) = Am(Pi) + A 2 n(P 2 ) + A 3 n(P 3 ) 

where \,i — 1,2,3 are barycentric coordinates. As pointed out above, n £ Nh if and only if 
|n| = 1 at all the mesh nodes. Thus we have 

\n(x)\ = |A 1 n(P 1 ) + A 2 n(P 2 ) + A 3 n(P 3 )| 

< A 1 |n(P 1 )| + A 2 |n(P 2 )| + A 3 |n(P 3 )| 
= Ai + A 2 + A3 
= 1 

If |^(a;)| = 0, then obviously the conclusion is true. In the following, we assume > 0. 

Let h = n(x)/\n(x)\, then \n\ = 1. We have 

(|F| 2 -(l-a)|F T n| 2 ) 
= (\F\ 2 -(l-a)\n(x)\ 2 \F T n\ 2 ) 

>(\F\ 2 -(l-a)\F T n\ 2 ) 
>a|P| 2 

where in the last step we have used Proposition |4j □ 
Remark Similar arguments work for 3D, as well. 



Theorem 12. There exists solution to the discrete minimization problem (51). 



Proof. Take any (uh, nh) £ A(uoh,noh), we have by Lemma 11 

n(tt fc ,n fc )> / a\(I + Vu h )\ 2 + b\Vn h \ 2 dx 
Jn 



~ \\f\\Li(n)\\Uh\\Li(n) - \\g\\Li(T 2 )\\uh\\u>{T 2 ) 

> 



/ a\(I + Vu h )\ 2 + b\Vn h \ 2 dx 

\\f\\h(n)+4u h \\ 2 L 2 {n) ) - [ j||fl||L2 ( r 2 ) +4u h \\h {r2 



> [ C 1 \(I + Vu h )\ 2 +C 2 \Vn h \ 2 dx-C 3 
Jn 



where e > is small and Ci > 0,i — 1,2,3 are constants and we have applied the general- 
ized Poincare inequality ([3], p281) and the Trace Theorem p258) in the last step. Thus 
IL(uh,nh) — ► 00 as or goes to 00. Therefore its minimum must be achieved at a 

bounded subset of A(uQh, ngh). 

On the other hand, since A is the intersection of reciprocal images of of the continuous 
functions gj, it's a closed set. 

Now we are minimizing a continuous function H(uh, nh) on a closed bounded finite- dimensional 
set, by the Weierstrass Theorem, we can find (uh, nh) £ A(uoh, n$h) minimizing II on A(uQh, noh). 

□ 

Remark This theorem is the discrete version of Theorem [7| Since we are dealing with finite 
dimensional spaces, the proof is much simpler. Plus, in the continuous case we were only able to 
prove Theorem [7] for 2D, but in the discrete case, it's easy to see that the proof above applies to 
3D, too. 
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3.2 Equilibrium equation and its linearization 

Similar to the continuous case, we start with the following energy functional with Lagrange 
multipliers 

£(u h ,n h ,p h ,X h )= [ (\F h \ 2 -(l~a)\Fln h \ 2 )+b\Vn h \ 2 
Jn 

- Ph (det(F h ) - 1) + \{ir h (n h ,n h ) - 1) (59) 
f u h - g u h da 



where F^ = I + Vit/j is the deformation gradient. 

Taking the first variation variation on ( 59 1 gives us the Euler-Lagrange equations (equilibrium 
equations) for the constrained minimization problem (57 1, and taking another variation on the 
Euler-Lagrange equations gives us the linearized problems. 

Thus the equilibrium equations are, find (u h ,n hl p h . \ h ) G W Mo |r xV Ml) | ri x V h x T4,A |ri: 
such that 



/ 2(F h : Vv-(l-a)(F£n h ,Vv T n h ))- Ph ^^ :Vv 
Jn u^h 

- I f-v- [ g-vda (60) 
Jn Jt 2 

[ -2(1 - a){F^n h ,F^m) + 2b(Vm,Vn h ) +2\ h n h (n h ,m) (61) 
Jn 

= / -q(6ebF h -l) (62) 
= / IMr h ((n h ,n h ) - 1) (63) 











for any test functions (v, m, q, fx) £ W /li0 |r x V h .. |ri x Vh x ^ft,o|ri i where Fh = I + Vw/i. 

The corresponding linearized problem is, for a given (it^, n^p/j, A/>) 6 W;, , Uo \r X „ | ri x 
Vh x V hM \ Tl , find (w,l,o,j) e W M |r x V M[ri x V h x V M | ri such that 

ai{w,v) + a 2 (l,v) + bi(o,v) = Li(v) (64) 

a 2 (m, w) + a 3 (l, m) + 6 2 (7, m) = L 2 (m) (65) 

6i(g,u7) = L 3 ( 9 ) (66) 

fcOM) = L 4 (M) (67) 

is true for any (v,m,q, n) G W/ t |r x Vft^iri x V/j x V/^oiiv Here ai, a 2 and b\ are the same as 
in the continuous case, while 

a 3 (l, m)= -2(1 - a)(F%l, F^m) + 26(Vm, VZ) + 2\ h n h (l, m) (68) 

and 

b 2 (^i,m) = / 2fnr h (n h ,m) (69) 
Jn 



3.3 Well-posedness of the linearized system 

As in the continuous case, the well-posedness of the linearized system can be reduced to the 
verification of the inf-sup conditions for the bilinear forms 6i(-, •), b 2 (-, •) and a(-, •). 

Since we used the Taylor-Hood element P 2 — P± for the (u,p) combination, the inf-sup con- 
dition for &i(-, •) is at least satisfied at the strain- free state (Proposition 6.1 of [3]). Since at the 
stress-free state the deformation gradient F is a constant, by change of variables, it's easy to 
verify that the inf-sup condition for b\(-, •) is satisfied at the stress-free state, as well. 
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For the inf-sup condition for &a(v)> we can use the results from [12]. A slight modification of 
the proof in [12] gives us the following result: 

Theorem 13. Assume n e H*,^ (fi) f| W 1 ' 00 ^), and n h e V h>no \ ri satisfies \n h \ > C > 

and \\rih — 7T/jn||i < 7/I log(/i)| 1 / 2 . TTien we can find a positive constant P2, independent of h, 
such that 

. (TTh[n h ■ m],fi) 

inf sup — — — — > p 2 (70) 
^V h , 0]Tl m£Vk, | ri ||/x||_i||m||i 

In general, analytical verification of the inf-sup conditions may not be easy. However, in the 
discrete case, we can relate the inf-sup values to the singular values of certain matrices. 
For the general inf-sup value in 

Ph = inf < sup b(q h ,Vh)> (71) 

q h £Ph,\\qh\\ = l {v h EV h ,\\v h \\=l J 

we have the following result ([3]). 



Theorem 14. The inf-sup value fih in (71) is equal to the smallest singular value of the matrix 



S 2 BT 2 ; where the matrices S, T, B are defined by the following equations 



dht^tfSq (72) 
v T Tv (73) 



v h 
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b{q h ,v h ) = q T Bv (74) 

and q, v are the degrees of freedom of qh and Vh respectively. 

To verify the inf-sup condition or ellipticity condition for a(-, •) on K.er(£&h), we are interested 
in calculating the inf-sup value (3h in 

$ h = inf sup a(u h ,v h ) on Ker(^ h ) C Y h (75) 

l|Mhll = l ||« h ||=l 

and the ellipticity constant 6th in 

a h = inf a{v h ,v h ) on Ker{@ h ) C Y h (76) 
ll"hll=i 

as well. 



Theorem 15. The inf-sup value fth in (75) is equal to the smallest singular value of the matrix 
A\, while the ellipticity constant &h in (76) is equal to the smallest eigenvalue of the matrix A\, 
where A\ is the lower right (n — m) x (n — m) corner of the matrix Q T T~ 1 / 2 AT~ 1 / 2 Q . Here n 
is the dimension ofYh, Tn is the dimension ofVh, and the matrices T, A, B are defined by the 
following equations 

\\v h \\ 2 = v T Tv 
a(u h ,v h ) = u T Av 
v h e Ker(^ h ) ^ Bv = 

where u, v are the degrees of freedom of Uh and Vh respectively. Here B is assumed to be full-rank 
and the matrix Q is defined by the QR decomposition of (BT~ 1 / 2 ) T 



{BT -l/2 )T = q 
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Proof. First, let x — u, y — v then we have 



k= mf. sup X- AV r^- (77) 

a;GKer(B) 3/eKer (B) V X 1 Xy'y 1 y 



where B = BT^ 1 / 2 , and I = T~ X I 2 AT~ X I 2 . 

Since i? is full-rank, the matrix R € M mxm in the QR decomposition 



U< - Q[f) (78) 



= ( ™* ) (79) 



is non-singular. Let 

«'—(, 

where w x e M m and z x e R n " m . Then it's easy to verify that 

a; e Ker(B) &w x = Q 
Thus there will be no constraint on z x . Therefore 

inf sup = infsup Z * Al ^ (80) 

xGKcr(B) xeK cr(S) VX T Xy/y T y z * z y z^Z x Jz^Z y 

where A\ is the lower right (n — m) x (n — to) corner of the matrix Q T AQ. Thus ^ is the smallest 
singular value of the matrix A\. 

Similar arguments reduce &h to the smallest eigenvalue of the matrix A\. □ 

Remark The matrix B is full-rank if and only if the operator S8h is onto, which is true if and 
only if the corresponding inf-sup condition holds. 

To compute the inf-sup value for 6 2 , we need to calculate the ifp 1 norm for any function Vh in 
^i,o|rv By the Ricsz Representation Theorem, we can find w € 7?g| ri , such that = 
The H 1 norm of v can be approximated by the H 1 norm of Vh, which is the L? projection of 
v e ^o|ri mto ^»,o|ri • Let {(fi, i = 1, n} be a basis of V/^oirj • We want to assemble the matrix 
S such that 

IKIIjy- 1 ~ \\vh\\m = vTSv 
r i 

where v e K" is the degree of freedom for t^. 

Theorem 16. The matrix S = AB~ X A, where the matrices A and B satisfy 

\\v h \\ L 2 = v T Av 
WviiWh 1 = v T Bv 

for any v h in Vf lt o\r 1 , and v <G K™ is i/ie degree of freedom for v h . 

Proof. Let / : ifp 1 — > V/j.o|ri be the map taking to u/,, and let fa — f(<fi). It's easy to see 
that 

Sij = {fa,fa)m (81) 

By definition of fa , we have 

J <Pi<Pj = J DfaDtpj + J <pnpj Vl<i,j<n (82) 
Since fa e Vf l fi\r 1 , we can write 

fa — ^ Gjfe^fc 
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Substituting it into ( 82 ) gives 

J (filfj = ^ G ik {^J Dip k ■ Dlfj + J 



Pklfij 



k 

That is A = GB. or G = AB^ 1 . Therefore 



S = ({<Pi,<f>j)Hi) 

= (Dtp^Dtpj) + (<pi,<fj) 

= ^ GipG]q [(D(fp, Dip q ) + (<p p , tp q )\ 



p,q 

y^^GipBpqG^j) 

GBG T 

(AB~ 1 )B(B~ 1 A) 
AB~ l A 



□ 



3.4 Existence and Uniqueness of the Lagrange multipliers 

Le Tallec proved in [15] that if the inf-sup condition is satisfied at it/, , then there exists a unique 
Ph such that (uh,Ph) is a solution of the discrete equilibrium equations of the incompressible 
elasticity. In this section, we will prove similar results for our elastomer problem. 

The proof uses the following result of Clarke [6] on constrained minimization problems, where 
dg(x) is the generalized gradient introduced in [5]. 

Theorem 17. Denote J a finite set of integers. We suppose given: E a Banach space, go, 
gj(j G J) locally Lipschitz functions from E to K, and C a closed subset of E. We consider the 
following problem: 

Minimize go(x) 

subject to x G C, 9j{x) — 0, Vj G J. (83) 



If x is a local solution of (83) then there exist real numbers ro, Sj not all zero, and a point £ in 



the dual space E' of E such that: 

£ G r Q dg (x) + J2 s jdg 3 (x), -f G N c (x), (84) 
j 

where N c (x) is the normal cone at C in x, dgj is the generalized gradient of gj(x). 

Theorem 18. Suppose (uh,nh) G Kh x N^, and at (uh,Tih), the inf-sup conditions for b\ and 
62 are both satisfied. Then there exist a unique ph G Vh and a unique A/ t G Vh.\ \r such that 



(uh,rih,Ph, ^h) is a solution of the discrete equilibrium equations \6(fy-(63). 
Proof. Denote 

E = C = (W M |r + u oh ) x (V hi0 | ri + ™oh) (85) 
g {x) = Il(v h ,m h ), gj(x) = gj(v h ,m h ) (86) 



where the functions <?j arc defined as in (55 1. It's easy to see that 



N c (x) = N E (x) = (0,0). (87) 
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Notice that 

dU(u h ,n h ) c {Dgl+Dgl} 
where Dgl and Dgl are in [(W fc)0 |r + u oh ) X (V M | ri + n oh)]* ■ We have 

Dgl(u h ,n h ) ■ (v h ,m h ) ■■ 



h(v h ) 
f2(m h ) 



where 



and 



Also, 



fi(v)= [ 2(F h :Vv-(l-a)(F^n h ,Vv T n h )) 
Jq 

him) = [ -2(1 - a)(F^n h ,F^m) + 2b(\7m,\7n h ) 
Jn 



Dgl(u h ,n h ) ■ (v h ,m h ) = 



Jnf- V h- Ir 2 9 ■ v h da 




We also note that gj(vh, mn) is continuously differentiable on (W^ |r + Uoh) x (V^^iri 
and that 



<% = {Dgj}, 
Dgj(u h ,n h ) ■ {v h ,m h ) 
tor 1 - ' i • X 
Dgj(u h ,n h ) ■ (v h ,m h ) = 



/n-Vi^(J + V«h):V« h 



J n ipj^ N Tr h [2n h ■ m h ]da 



for N + 1 < j < N + M 

Therefore applying Theorem |17[ we have: 

There exists real numbers ro, Sj not all zero, such that 

N+M 

er dU(uh,nh)+ ^2 Sjdgj(uh,n h ) 



Using ( 88 ) and ( 89 1 , equation ( 92 ) becomes 



N+M 



r {Dgl(u h ,n h ) + Dg%(u h ,n h )} + ^ SjDgj(u h ,n h ) = 
in [(W M |r + u oh ) x (V /lj0 |r 1 + n oh )]* . 



Suppose now tq = 0. By linearity, and using equations (90), (91), we can then transform 
get 



ddet 
dF 



(I + Vu h ) : Vv h dx = 0, Vv h G W M |r 



/ y2 s N+j^P] n h {2n h ■ m h ]dx = 0, Vm ft <E V^ri) 
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Since at least one Sj is nonzero, at least one of the equations (94 1 or (|95|) is in contradiction with 



the inf-sup conditions. Thus r$ cannot be zero. We can then divide ( |93[ ) by tq to get 

M 

Dgl(u h ,n h ) ■ {v h ,m h ) + l/r ^ SjDgj(uh, nn) = -Dgl(u h ,n h ) ■ (v h ,m h ), 



That is 



V(v h> m h ) e [(W M |r +u oh ) x (V M | ri + n oh )] 
ddet 



A(v h ) 



f d dot P f 

I Ph^Eri 1 + Vm >0 ■ Vv h dx = / / • v h + / g v h da 
Jn <Jt j n j r 



h{ rn h)+ \ X h TTh[2n h ■ m h ]dx = 
Jn 



where we have denoted 



JV 



Ph 



u 



^SN+jipj ] /r Q 



(96) 

(97) 
(98) 

(99) 
(100) 



conclude that (uh,rth,Ph, ^h) is a solution of (60)-(63). 



Equations (97l and (98) are precisely ( 60 )-( 61 1. Since we also have (uhj-n^) e K] L x N^, we 



Finally, if there were two different p^, their difference would violate the inf-sup condition for 
bi; and if there were two different A^, their difference would violate the inf-sup condition for &2- 
So we have the uniqueness of both ph and A/,. □ 



3.5 Implementation using FEniCS 

All the computations in this project were done using FEniCS ([ 
software for automated solution of differential equations. 



which is an open source 



We want to solve the equilibrium equations (60)-(63). Assume (it, n,p, A) are the degrees 
of freedom for (u^, n^Ph, A/j). Replacing the test functions (v,m,q, fj.) by the corresponding 



basis functions, the equilibrium equations (60)-(63) can be looked as a nonlinear equation for 
(w,n,p, A): 

&(u,n,p,\) = (101) 
Since it's nonlinear, we can use Newton's method to solve. The iterations are, for integers k > 0, 
D,^(u k ,n k ,p k ,X k ) 



D{u,n,p, A) 
It can be verified that the matrix 



-(Au k , An K ,Ap k ,AX k ) = -&(u k ,n k ,p k , \ k ) 



A 



D,^(u k 1 n k ,p k ,\ k ) 
D(u,n,p,\) 



(102) 



(103) 



is exactly the matrix corresponding to the left side of the linearized system ( 64 1-(|67[) (i.e. 



by replacing (w,l,o, 7) and (v,m,q, fx) by the corresponding basis functions, and replacing 
(uh,nh,Phi Xh) by the solution at step k). Each step of (102) corresponds to a linear variational 
problem, and can be solved by FEniCS . Thus the whole problem can be solved. 

There is a minor issue, though, about the procedure above. In both the equilibrium equations 



(60)-(63) and the linearized equations (64)-(67), there are terms with the interpolation operator 



7T/j, which is not supported in the FEniCS form language. Thus we can't directly input those 
bilinear and linear forms in the form file for FEniCS. We solved this issue by first letting FEniCS 
assemble the matrix A and the right-side vector b without the 7r^ terms, and then assemble the 
7T/J terms ourselves and update the matrix A and the vector b accordingly. It turns out that the 
7r/j terms are not very difficult to assemble. 
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4 Numerical Results 



The simulation setup is as in Figure [4] The reference domain is [0, AR] x [0,1]. The clamps 
constraint are imposed at X — and X = AR. At the beginning, the elastomer is pre-stretched 



Figure 4: The elastomer is clamped and pulled on both sides. 

to the stress-free state, and it's assumed that at this state the directors are aligned uniformly in 
the direction 

n = (0,l) T (104) 

The values of the displacement, pressure etc. at the stress-free state were mentioned in section 
|2.3| For completeness, we rewrite them here. The displacement at the stress-free state is given 

by 

u x = (a 1/4 — 1)(X — 0.5AR) (105) 
u Y = (a~ 1/4 -l)(y-0.5) (106) 

The pressure at the stress-free state is given by 

p = 2a 1/2 (107) 

And A at the stress-free state is given by 

A = (l-a)/Va (108) 
Suppose we want the aspect ratio at the stress- free state to be AR„, we need to set 

AR = -^AR n (109) 
V a 

Notice that the whole problem is symmetric about both the vertical middle line X — 0.5AR 
and the horizontal middle line Y = 0.5. Hence we only need to do the computation at the upper- 
right quarter [0.5AR, AR] x [0.5,1]. By symmetry, we have the following Dirichlet boundary 
conditions: 

u x = on X = 0.5AR (110) 
u Y =0 on 7 = 0.5 (111) 
n = (0, 1) T on X = 0.5AR and Y = 0.5 (112) 

Also, we model the "clamped pulling" by specifying the following Dirichlet boundary conditions 
at X = AR: 

u x = 0.5AR[a 1/4 (l + Mi) - 1] (113) 
u Y = (a- 1/4 -l)(y-0.5) (114) 
n=(0, 1) T (115) 
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where t G [0, 1], and 1 + M is the largest elongation factor (ratio of the length of the elastomer 
after stretching and the length of the elastomer before the stretching). We slowly increases t from 
to 1, to increase the numerical stability of the Newton's method. To make sure that \n\ = 1 at 
all the mesh nodes, we need to impose Dirichlet boundary conditions for A at the same boundary 
as n. Therefore we impose 

X=(l-a)/y/a on X = 0.5AR, X = AR and Y = 0.5. (116) 

We took a = 0.6, b = 0.0015, AR„ = 1, M = 0.4, with "time" step At = 0.01. We have used 
uniform triangular meshes in the simulations. For example, Figure [5] shows the mesh with mesh 
size h = 2~ 5 (notice that the mesh is 16 x 16, but recall that this computational domain is only 
1/4 of the true domain) . 




ure 5: The uniform meshes we have used. What's shown here is the mesh with mesh size h = 

Table [l] summarizes the numerical errors. We can see that the L 2 errors of Uh and rih, and 
H~ x error of Ah are relatively small, while the H l errors of U/, and rih are relatively large. Table 



h 


2~' 2 


2- :i 


2 _4 


2" 5 


\\ u h 


- Uh/2 





3.49E-03 


1.91E-03 


8.39E-04 


2.69E-04 


\\ u h 


~ U h/2 


1 


5.14E-02 


3.77E-02 


2.02E-02 


7.66E-03 


\\ n h 


- n h / 2 





2.32E-01 


9.70E-02 


3.05E-02 


8.25E-03 


\\ n h 


- nh/2 


1 


2.31E+00 


1.91E+00 


1.19E+00 


6.23E-01 


\\Ph 


- Ph/2 1 





1.68E-01 


7.93E-02 


2.38E-02 


8.99E-03 


IIAh- 


~ ^h/2\\ 


-1 


1.15E-02 
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1.51E-03 


5.22E-04 



Table 1: The numerical errors. 



[2] summarizes the convergence rates. We can see that the convergence rates are relatively fast, 
except the H l errors of Ut and rih- 

Next, let's check the results of the inf-sup tests. Table[3]lists the inf-sup values at t = 0, where 
s(^/|Ker(^)) is the inf-sup value of ^|Ker(^), while e{j^\Kej{3§)) is the ellipticity constant of 
si/\Ker(S§). We can see that the inf-sup values of b\ and hi don't change much with the mesh. 
This implies that the inf-sup conditions for b\ and 62 probably are always satisfied regardless of 
the mesh. Next, we can see that the inf-sup values for s(^|Ker(^)) are positive for all the 4 
meshes, which means the inf-sup conditions are satisfied thus the discrete saddle point systems 
are well-posed for all the 4 meshes. However, it seems that the inf-sup values s(jz/|Ker(^)) tend 
to as h goes to 0. Finally, we can see that the ellipticity constants e(si/\Ker(£§)) are negative for 
the given 4 meshes. This means the ellipticity conditions are not satisfied at the initial stress-free 
state. This is consistent with our analytical observation. Table [4] lists the inf-sup values at the 



21 



h 


2" 3 


2 _4 


2" 5 


log 2 (\\u 2 h 


- u h 


lo/ 


\u h 


~ u h/2 


lo) 


0.87 


1.18 


1.64 


log 2 (||«2/i 


- u h 


li/ 


\u h 


~ u h/2 


ll) 


0.45 


0.90 


1.40 


log 2 (||n 2 h 


- n h 


lo/ 


\n h 


- n h / 2 


lo) 


1.26 


1.67 


1.88 


log 2 (\\n 2 h 


- n h 


ll/ 


\n h 


- n h/2 


ll) 


0.27 


0.69 


0.93 


log 2 dlm 


-Ph 


lo/ 


\Ph 


~Ph/2\\ 


3) 


1.08 


1.74 


1.41 


log 2 (||A 2 /i - 


- A fc ||_ 


-1/ 


\Xh 


- \/2\ 


-l) 


1.38 


1.55 


1.54 



Table 2: The convergence rates. 
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0.5836 


0.5875 
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2.0000 


2.0000 


2.0000 


2.0000 




3.60E-03 


2.70E-04 


5.69E-05 


1.62E-04 


e(^Ker(^)) 


-3.60E-03 


-1.27E-02 


-1.78E-02 


-1.21E-02 



Table 3: The inf-sup values at t = 0. 



final state t — 1 (elongation factor s = 1.4). We can see that the inf-sup values for b\ and 62 still 
don't vary much with respect to the mesh size h. This suggests that the inf-sup conditions for b\ 
and 6 2 are probably satisfied for all the u h and n h solutions during the pulling process. On the 
other hand, we still see that s(£/\Ker(&)) are positive, while e(&/\Kei(33)) are negative, which 
means the inf-sup condition for s(«e/|Ker(J?)) are satisfied, while the ellipticity condition are not. 
Similar to t = 0, although the inf-sup values s(^\Ker(£§)) are positive for the 4 given meshes, 
they do seem to converge to zero as the mesh size h goes to 0. 
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0.6549 


0.6431 
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b 2 


1.9967 


1.9503 


1.9065 
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1.20E-03 


5.82E-04 


4.88E-05 


e(j*lKer(#)) 


-2.91E-03 


-2.58E-03 


-5.82E-04 


-4.88E-05 



Table 4: The inf-sup values at t = 1. 

Figure [6] shows the graph of nominal stress vs. strain, where we have taken the stress-free 
state as the reference configuration. We can clearly see a plateau when the strain is in (0.10, 0.22). 
Comparing with Figure [2j we can say that we've recovered the semi-soft elasticity phenomenon. 
The plateau regions in Figure [2] are more flat than ours, because those graphs are for long and 
thin elastomers, while our graph is for the elastomer with the aspect ratio AR„ = 1. 

Figure [7p] and [9| show the directors configuration at the start, bottom, and end of the plateau 
region (elongation factor is about 1.10, 1.17 and 1.22 respectively). We can see that the directors 
just start to rotate at the start of the plateau region (Figure [7 1 , while many of the directors have 
finish the rotations at the end of the plateau region (Figure 9l). Also, we can see from Figure 
7p and |9] that during the plateau interval (0.10,0.22), the elastomer domain is mostly blue (low 
BTW energy); while after the plateau region (Figure [l0| , the red part (hight BTW energy) starts 
to dominate the elastomer domain. This means by rotating the directors, the elastomer maintains 
relatively low BTW energy during the elongation; while after most of the directors have already 
finished the rotation, the BTW energy starts to increase. This suggests that the plateau region 
in the stress-strain graph is probably caused by the rotation of the directors. This agrees with 
the theory of liquid crystal elastomers [TH] • 

Figure 10 shows the final configuration of the directors (elongation factor s — 1.4). We can 
see that most of the directors have finished the rotation in the final state. 
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Figure 6: Nominal stress vs. strain. The plateau region corresponds to strains in the interval 
(0.10,0.22). 
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Figure 8: The directors at the bottom of the plateau region, elongation factor s = 1.17. 




Figure 9: The directors at the end of the plateau region, elongation factor s = 1.22. 



24 




Figure 10: The directors at the final state, elongation factor s = 1.4. 



5 Conclusion and Discussion 

We have designed a numerical algorithm to simulate the 2D BTW+Oseen-Frank model of liquid 
crystal elastomer. We aimed to recover numerically the experimentally observed phenomena such 
as the stripe-domain phenomenon, the semi-soft elasticity etc. The existence, well-posedness, 
uniqueness and convergence etc. of the numerical algorithm have been investigated. 

We have successfully captured the semi-soft elasticity of liquid crystal elastomer, and found 
it to be related with the rotation of the directors. However, we didn't see the stripe-domain 
phenomenon. We believe this is mainly because the b value (the coefficient of the Oseen-Frank 
energy) 0.0015 is relatively large comparing to true value of b in practice, which we estimate 
to be probably less than 10~ 6 . Since the b value is relatively large, the Oseen-Frank energy 
dominates the BTW energy, while the latter of which is the main driving force for the occurrence 
of stripe domains. On the other hand, to observe the stripe domain phenomenon, we should 
use meshes fine enough to resolve the stripes. For example, in the experiment of Zubarev and 
Finkelmann ([20]), the width of stripe domain divided by the width of the elastomer is around 
15/j,m/5mm = 0.003, which is much smaller than the mesh size 2 -5 of our finest mesh. However, 
very small b values require small "time" steps to be stable, and this together with very fine meshes 
would be computationally very expensive. We have tried 512 x 512 uniform mesh, and the program 
ran out of memory. In the future, when the machines have much greater computing power, we 
can try very small b values and very fine mesh, and see whether the stripe domain occurs. Instead 
of using uniformly fine mesh, we can also use adaptive mesh to reduce the computational cost. 

Another direction worth trying is to use Ericksen or Landau-de Gennes' model instead of the 
Oseen-Frank energy in the numerical simulation. Oseen-Frank energy only allows point defects, 
while Ericksen or Landau-de Gennes' model allows line and surface defects, as well ([IS]). The 
stripe domains might have line or surface defects in the transition area between the stripes, thus 
using Ericksen or Landau-de Gennes' model might have a better chance of capturing the stripe 
domain phenomenon. 
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